arXiv:1504.06015vl [cs.IT] 23 Apr 2015 


Super-Resolution of Mutually Interfering Signals 


Yuanxin Li 

Department of Electrical and Computer Engineering 
The Ohio State University 
Columbus, Ohio, 43210 
Email: li.3822@osu.edu 


Yuejie Chi 

Department of Electrical and Computer Engineering 
The Ohio State University 
Columbus, Ohio, 43210 
Email: chi.97@osu.edu 


Abstract —We consider simultaneously identifying the mem¬ 
bership and locations of point sources that are convolved with 
different low-pass point spread functions, from the observation 
of their superpositions. This problem arises in three-dimensional 
super-resolution single-molecule imaging, neural spike sorting, 
multi-user channel identification, among others. We propose a 
novel algorithm, based on convex programming, and establish 
its near-optimal performance guarantee for exact recovery by 
exploiting the sparsity of the point source model as well as inco¬ 
herence between the point spread functions. Numerical examples 
are provided to demonstrate the effectiveness of the proposed 
approach. 

Index Terms —super-resolution, parameter estimation, atomic 
norm minimization, mixture models 


I. Introduction 

In many emerging applications in applied science and 
engineering, the acquired signal at the sensor can be regarded 
as a superposition of returns from multiple channels (or 
users), where the return from each channel is governed by 
the underlying physical held that produced it, e.g. the Green’s 
function, the point spread function, the signature waveform, 
etc. The goal is to invert for the held parameters of each 
channel that produced the acquired signal which rehects the 
ensemble behavior of all channels. 

Mathematically, consider the acquired signal, y{t), given as 

I 

Xi{t)*gi(t) = ^ 

Z —1 2 = 1 

where * denotes convolution, Xi{t) = — Uk) is 

the point source signal observed through the ith channel, gi (t) 
is the point spread function of the ith channel, respectively, and 
I denotes the total number of channels. Eor the ith channel, 
let Tik € [0,1) and aik G C be the location and the amplitude 
of the fcth point source, 1 < fc < Ki, respectively. In typical 
applications, we are interested in resolving the point sources 
of each channel at a resolution much higher than that of the 
acquired signal y{t), determined by the Rayleigh limit, or in 
other words, the bandwidth of the point spread functions. 

The proposed model dTli occurs in a wide range of practical 
problems, ranging from three-dimensional super-resolution 
single-molecule imaging, to spike sorting in neural recording 
and DNA sequencing, to multi-user multi-path channel iden- 
tihcation in communication systems, and many others. 


yit) = X 


tkikgi{t Tik) 


( 1 ) 


Three-dimensional super-resolution imaging: By employing 
photoswitchable fluorescent molecules, the imaging process 
of stochastic optical reconstruction microscopy (STORM) H] 
is divided into many frames, where in each frame, a sparse 
number of fluorophores (point sources) are randomly activated, 
localized at a resolution below the diffraction limit, and 
deactivated. The final image is thus obtained by superimposing 
the localization outcomes of all the frames. This principle 
can be extended to reconstruct a 3-D object from 2-D image 
frames 0, by introducing a cylindrical lens to modulate the 
ellipticity of the point spread function based on the depth of the 
fluorescent object. Therefore, the acquired image in each frame 
can be regarded as a superposition of returns from multiple 
depth layers, where the return from each layer corresponds 
to the convolution outcome of the point sources in that depth 
layer with the depth-dependent PSE, as modeled in ([TJ. The 
goal is thus to recover the locations and depth membership 
of each point source given the image frame that records the 
returns from all depth layers. 


Spike sorting for neural recording: Neurons in the brain 
communicate by bring action potentials, i.e. spikes, and it 
is possible to hear their communications through a single or 
multiple microelectrodes, which record simultaneously activ¬ 
ities of multiple neurons within a local neighborhood. Spike 
sorting El, thus, refers to the grouping of spikes according 
to each neuron, from the recording of the microelectrodes. 
Interestingly, it is possible to model the spike fired by each 
neuron with a characteristic shape i). The neural recording 
can thus be modeled as a superposition of returns from 
multiple neurons, as in ([TJ, where the return of each neuron 
corresponds to the convolution of its characteristic spike shape 
with the sequence of its bring times. A similar formulation also 
arises in DNA sequencing 0. 

Multi-path identification in random-access channels: In 

multi-user multiple access model 0, the base station receives 
a superposition of returns from active users, as in dD, where 
the received signal component for each active user corresponds 
to the convolution of its signature waveform with the unknown 
sparse multi-path channel from the user to the base station. 
The goal is to identify the set of active users, as well as their 
channel states, from the received signal at the base station. 


A. Related Work and Contributions 


algorithm is demonstrated in numerical experiments. 


There’s extensive research literature on inverting ([T]i when 
there is only a single channel with / = 1, where conven¬ 
tional methods such as matched filtering, MUSIC El, matrix 
pencil 0, to more recent approaches based on total variation 
minimization El, can be applied. However, these approaches 
can not be applied directly when multiple channels exist in 
the observed signal, due to the mutual interference between 
channels. To the best of the authors’ knowledge, methods for 
inverting ([T]) with multiple channels have been extremely lim¬ 
ited. In 0, Eoi, im, sparse recovery algorithms have been 
proposed to invert ([T]i by assuming the source locations lie on 
a fine-grain grid, which is mismatched from the actual source 
locations, resulting in possibly performance degradation ifT^ . 
Even when all the point sources indeed lie on the grid, existing 
work m shows that the sample complexity needs to grow 
logarithmically with the size of the discretized grid, which is 
undesirable. The continuous basis pursuit algorithm 03 can 
be applied to retrieve source locations off the grid, but lacks 
performance guarantees. 

In this paper, we study the problem of super-resolving O 
when there’re two channels, i.e. 1 = 2. We start by recognizing 
that in the frequency domain, the observed signal can be re¬ 
garded as a linear combination of two spectrally-sparse signals, 
each composed of a small number of complex sinusoids. We 
then separate and recover the two signals by motivating their 
spectral structures using atomic norm minimization, which 
has been recently shown as an efficient convex optimization 
framework to motivate parsimonious structures m, El, 03, 
as well as satisfying the observation constraints. 

The separation and identification of the two channels, using 
the proposed algorithm, denoted by convex demixing, is made 
possible with two additional conditions. The first condition 
is that the point source models satisfy a mild separation 
condition, such that the locations of the point sources are 
separated by at least four times the Rayleigh limit; this is 
in line with the separation condition required by Candes 
and Fernandez-Granda El even with a single channel. The 
second condition is that the point spread functions of different 
channels has to be sufficiently incoherent for separation, 
which is supplied in our theoretical analysis by assuming 
they’re randomly generated from a uniform distribution on the 
complex circle. We demonstrate that, as soon as the number 
of measurements, or alternatively, the bandwidth of the point 
spread functions, is on the order max(A"i, ^^ 2 ) log(iTi + K 2 ) 
up to logarithmic factors, the proposed algorithm, denoted by 
convex demixing, recovers the locations of the point sources 
for each channel exactly, with high probability. Since at 
least an order of Ki + K 2 measurements is necessary, our 
sample complexity is near-optimal up to logarithmic factors. 
Moreover, the point sources can be recovered from the dual 
solution of the proposed algorithm, without estimating or 
knowing the model order a priori. Our proof is based on 
constructing two related polynomials that certify the optimality 
of the proposed algorithm. The effectiveness of the proposed 


B. Organization 

The rest of this paper is organized as follows. The proposed 
algorithm based on convex programming and its performance 
guarantee are provided in Section HI] The proof of the per¬ 
formance guarantee is sketched in Section |III| Numerical 
experiments are shown in Section |IV] and we conclude the 
paper in Section |V] Throughout the paper, (•)^ and (•)* 
denote the transpose and Hermitian transpose respectively, (•) 
denotes the element-wise conjugate, Tr (•) denotes the trace 
of a matrix, and toep {u) denotes a Toeplitz Hermitian matrix 
with u as its first column. 


II. Super-resolution oe Mutually Interfering 
Signals 


Denote the discrete-time Fourier transform (DTFT) of gi{t) 


9 i,n — 




( 2 ) 


which satisfies gi^n = 0 whenever n ^ = 

{—2M,..., 0,..., 2M}, where 2M is the cut-off frequency of 
band-limited gi (t). Typically, M is determined by the physics, 
such as the aperture of the imaging device or the steering array. 
Taking the DTFT of ([T]i with J = 2, we obtain 


yn = I 

' —00 

■K, \ 

^-j27rnTik 


— 9^,n ' ( ^ ^ 
\fe=l 


' K 2 


■g2,n- ^a2fce-^'2™ 


V 




(3) 


with n G fljvr- The measurements yn’s in (|3 can be considered 
as a linear combination of two spectrally-sparse signals, with 
gin’s determining the combination coefficients. Multiplying 
both sides of (|3 with g^J^, and with slight abuse of notation, 
we rewrite 0 into a vector form: 


y = + gOx2, (4) 

where y = [ 1 /- 2 M, • ■ • , 2/o, • ■ •, 2/2 m]’" S C^^+\ g = 
[ 5 - 2 M, • ■ ■ ,50, ■ ■ ■ ,ff2M] G with = P2,n/3l,n, 

and © denotes the Hadamard element-wise product operator. 
Furthermore, let x\ = [a:i,_ 2 M, • ■ •, • ■ •, S 

(^4M+l ^ [X2 -2M, ■ • ■ , 2 : 2 , 0 , ■ ■ ■ , X2,2mY' ^ 0^*^+^ 

denote two spectrally-sparse signals, each composed of a small 
number of distinct complex harmonics, represented as 

Ki 

®2 = X] (r2fc) € 

fc=l 

where Ki and K 2 are the spectral sparsity levels of each 
signal. The atom c (r) is defined as 

c(t) = . 1 . g-j2TTi2M)T ^ 




which corresponds to a point source at the location r € [0,1]. 
Further denote the sets of point sources in x* and a:^ by Ti = 
{tii, ..., Tiif J and T 2 = {r 2 i,..., T 2 k. 2 } respectively. The 
goal is thus to recover Ti and T 2 , and their corresponding 
amplitudes, from (IHi. 

Define the atomic norm iflTl . ifTSl of a; G c4M+i 
respect to the atoms c(r) as 

ll=®IU = inf , i y] lofel I at = Vafec(rfc) I , 

afceC,TfcG[ 0 ,l) ^ J 

which can be regarded as the tightest convex relaxation of 
counting the smallest number of atoms that is needed to 
represent a signal x. Therefore, we seek to recover the signals 
ati and X 2 by motivating their spectral sparsity via minimizing 
the atomic norm, with respect to the observation constraint: 

{xi,X 2 } = argmin||a;i|U + ||a; 2 |U, 

^l,X2 (6) 

s.t. y = Xi+gQx 2 . 


The above algorithm is referred to as convex demixing. Inter¬ 
estingly, (|6]l can be equivalently rewritten with the following 
semidefinite programming characterization Ea, which can be 
solved efficiently using off-the-shelf solvers, as 


Tr(toep(Mi))-f Tr(toep(M2)) , 

mm -:-h li -f 19 , 

Xi,X2.Ui,U2,ti,t2 AM + 1 


s.t. 


toep(Mi) xi 

>- 0 , 

toep (M2) X 2 

xl h 

— 5 

X*2 t2_ 


y = Xi+gQx2. 


>1 0, 


(7) 


Dehne the separation condition of point sources of each 
channel as 

Ai = min|Tifc - Tij\ , 

which is the wrapped-around distance on [0,1], and the 
minimum separation of all channels as A = min^ A^. Our 
main theorem is stated below. 


Theorem 1. Let M > 4. Assume that Qn = ’s are 

i.i.d. randomly generated from a uniform distribution on the 
complex unit circle with fn ^ W[0,1], and that the signs of 
the coefficients aik’s are i.i.d. generated from a symmetric 
distribution on the complex unit circle. Provided that the 
separation A > 1/M, there exists a numerical constant C 
such that 


functions, is on the order M = 0(max(iTi,/T 2 ) log(iTi -f 
K 2 ) log M), the proposed convex demixing algorithm recovers 
the locations of the point sources for each channel exactly, 
with high probability. This suggests that the performance of 
the convex demixing algorithm is near optimal in terms of the 
sample complexity. 

Remark 1. The point source separation condition A > 1/M 
is used as a sufficient condition in Theorem [T| to guarantee 
accurate signal demixing. It is implied in ||9l that a reasonable 
separation is also necessary to guarantee stable superres¬ 
olution. Interestingly, no separation between point sources 
is necessary across channels, as long as their point spread 
functions are sufficiently incoherent. 

Remark 2. Theorem [T] assumes g„’s are generated with 
uniformly random phase, which may be reasonable when 
5 „’s can be designed, such as the spreading sequences in 
multi-user communications. This assumption may be further 
relaxed. The proof procedure reveals that same results can be 
obtained as long as g„’s satisfy E [gn\ = E [ 9 ^^] = 0 and 
Cl < \gn\ < C' 2 . 

Remark 3. Both sign(aife) and sign(a 2 fc) are required to be 
randomly generated, which we believe are technical require¬ 
ments of the proof, and may be removed with hner proof 
techniques. 

Define the inner product of two vectors as {p, x) = x*p and 
the real-valued inner product as (p, a;)R = Re {x*p). Then the 
dual norm of H-H^ can be represented as 

IIpII^ = sup { p , x ) r = sup 

||a!:||^<l i-eiO.l) 

We have the dual problem of (| 6 ]l written as 

p = argmax (p, p)r, 

p (8) 

s-t. ||p||]4 < 1, WgQpfj^ < 1 , 

which comes from standard Lagrangian calculations. Construct 
two dual polynomials using the dual solution p: 

2M 2M 

Pir)= Q{r)= y Pn9ueP^^\ 

n——2M n——2M 

then the point source locations Ti and T 2 can be recovered 
without model order estimation, by identifying the parameter 
r such at |P(t)| = 1, and |Q(t)| = 1, respectively. 


2M 

y PneP' 

n——2M 


M > C max log' 
max{ifi,Ar2}log 


2 (M{Ki+K2) 


9 

K 1 +K 2 

V 


log 


M{Ki+K2) 


is sufficient to guarantee that x\ and ajJ are the unique 
solutions of (|6j) with probability at least 1 — y. 


Theorem [T] indicates that as soon as the number of mea¬ 
surements, or alternatively, the bandwidth of the point spread 


III. Proof Sketch OF Main Result 

In this section, we proceed to sketch the proof of Theorem 
[T| We first provide the optimality conditions using dual poly¬ 
nomials to certify the optimality of the solution of (|6ll. Illumi¬ 
nated by 121, na, where the dual polynomial is constructed 
by the squared Fejer’s kernel, we propose a construction of 
dual polynomials which are composed of a deterministic term 
and a random perturbation termed introduced by interference 
between channels. Finally, the remainder is to show that the 
constructed dual polynomials satisfy the optimality conditions 












with high probability when the number of measurements M 
is large enough. 


and 


A. Optimality Conditions using Dual Polynomials 

We can certify the optimality of the primal problem (|6]l 
using the following proposition. 


Proposition 1. The solution of ^ xi = and X 2 = X 2 
the the unique optimizer if there exists a vector p such that 
the dual polynomials P{t) and Q{t) constructed from it 


P(r) 

satisfy 


2M 2M 

^ Q (r) = ^ (9) 

n— — 2M n— — 2M 


{ P (rifc) = sign (aifc), Vrifc e Ti 
|P(t)|< 1, Vr^Ti 

Q (T 2 k) = sign (a 2 fc), VT 2 fc € T 2 ’ 
|Q(r)|<l, Vr^T 2 


( 10 ) 


where the sign should be understood as the complex sign. 


Proposition [T] suggests that P (r) and Q (r) are dual cer¬ 
tificates to show {x\^xf] is the unique primal optimizer. 
Therefore if we can find a vector p to construct two dual poly¬ 
nomials P{t) and Qir) in (|9]l satisfying (fTOl i. the proposed 
algorithm is guaranteed to recover the ground truth. 


B. Dual Certificate Construction 

Our construction of dual polynomials is inspired by 13, 
ina, based on use of the squared Fejer’s kernel. However, 
since the two dual polynomials are coupled together, the 
construction is more involved. Define the squared Fejer’s 
kernel 0 as 


Ki 

Q(t) = OiikKg (t - Tik) + 

t kT ( 14 ) 

+ - 'T2k) + ^ l32kK' (t - T2k) , 

fc=l fc=l 


where Tik € Ti and T 2 k S T 2 . It is straightforward to 
validate that there exists some corresponding vector p such 
that (fT3]) and (fT4l) can be equivalently written as (|9]). Set the 
coefficients ai = [an,..., (3i = [fin,..., , 

(y .2 = [a 2 i, - ■ ■ ,^ 2 X 2 ]^ and f32 = [fi 2 i, ■ ■ ■, fi 2 K 2 f by 
solving the following equations 


(Tifc) = sign (aifc), rifeSTi, 

P'iTlk)=0, TlfcGTl, 

Q {T 2 k) = sign (a 2 fe), T 2 k G T 2 , 

Q'{j2k)=Q, 'r2fc G T 2 . 


The above setting, if exists, immediately satisfies the first and 
third conditions in (fTOl i. The rest of the proof is then to, 
under the condition of Theorem [T] guarantee that a solution 
of (fTSl l exists with high probability, and that when existing, 
they satisfy the second and forth conditions in ( fTOb with high 
probability, therefore completing the proof. 

C. Validation of Constructed Dual Polynomials 
First we want to show that the solution of (fTSl) exists with 
high probability. Since the deterministic terms in constructed 
dual polynomials are well-conditioned if the point source 
separation A satisfies appropriate condition ifTSll . by showing 
the random perturbations are small, we can have the following 
proposition to guarantee the invertibility of (fTSl l with high 
probability when the number of measurements is large enough. 


2M 


n=-2M 


Sn^' 


j2'KnT 


( 11 ) 


where s _ J_ _LI) 

wnere Z^i=niax(n-M.-M) (^ I M I) (-^ I M m\)- 

The value of K (r) is nonnegative, attaining the peak at 

T = 0 and decaying to zero rapidly with the increase of the 

absolute value of r. 

We define Kg (r) and Kg (t) respectively as 

2M 

^n9n 

n^-2M 
2M 


ZM 

^9ir) = jT E “ " 


( 12 ) 


K-n 




^j27rnT 


-i=-2M 


We then construct the dual polynomials P (t) and Q (r) as 

ifi ifi 

(t) = E ('^ - Dk) + fiikK' (t - rife) 

fc=i fe=i ,,,, 

K 2 K 2 ^ ’ 

+ ^ a 2 kKg (r - T 2 k) + ^ P 2 kKg {r - T 2 k ), 


Proposition 2. Assume M > 4. Let 6 G (0,0.6376) and 77 G 
(0,1), then ( 1151 ) is invertible with probability at least 1 — p 
provided that 

M> ^max{7Fi,iT2}log(^E^^i^±iE^ . (I 6 ) 

Once the constructed dual polynomials are fixed with 
coefficients {ai,/3i, a 2 ,/32} from the solution of (fTsT l. the 
rest is to verify that |F’(t)| < 1, Vr ^ Ti and similarly, 
|(5(r)| < l,Vr (ji T 2 . Since the expressions for P{t) and 
Q(t) are very similar, it is sufficient to only establish the 
above for P{t). 

Denote 

^ (t) = E [ ^ (^)j , 

s/wm vvwm \ 

which is well-bounded M, 0, where the expectation is with 
respect to g and (r) is the Zth derivative of P (r). 

While the random perturbation in constructed P^^'> (r) in¬ 
troduced by interference effect can be verified small enough 
with high probability, we have the following proposition to 
uniformly bound the distance between (r) and (r). 









Proposition 3. Suppose A>^. If there exists a numerical 
constant C such that 


M > max | log^ 

UYSiK{Ki,K2}\og 
then 


2 fM{K^ + K2) 


€71 


K, +K2 


log 


M {Ki + K2) 


€7] 


1 


VW^\ 


.PW (r) - 




.pW (r) 


< e, 


M > C max log' 

max{pri,pr2}log 

then 


2 fM{K,+K2) 


V 

Ki + K 2 
1 


log 


M {Ki + K 2 ) 


V 


l^('^)l < 1: forr e [0, l]\Ti, 
with probability at least 1 — rj. 


\/t S [0,1), / = 0,1, 2, 3, holds with high probability at least 
l-rj. 

Then applying Bernstein’s polynomial inequality M and 
similar techniques in ifTSl Lemma 4.13 and 4.14], we have the 
following proposition. 

Proposition 4. Suppose A > ■^. If there exists a numerical 
constant C such that 



(a) M = 8 (b) M = 16 

Fig. 1: Successful rates of the convex demixing algorithm 
when (a) M = 8 and (b) M = 16. 


V. Conclusions 

We propose a convex optimization method based on atomic 
norm minimization to super-resolve two point source models 
from the measurements of their superposition, with each 
convolved with a different low-pass point spread function. It 
is demonstrated, with high probability, that the point source 
locations of each channel can be simultaneously determined 
perfectly, from an order-wise near-optimal number of measure¬ 
ments, under mild conditions. The proposed framework can be 
extended to handle more than two channels, whose details will 
be discussed elsewhere. 
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The proof of Theorem [T] is then complete since we have 
established that P{t) and Q{t) constructed in ( fOl ) and 
dl are indeed valid dual certificates under the condition of 
Theorem [T] 

IV. Numerical Experiments 

We carry out numerical experiments to validate the per¬ 
formance of the proposed convex demixing algorithm. For 
a fixed M, we vary the spectral sparsity level of the two 
channels as Ki and K 2 . For each pair of {Ki,K 2 ), we first 
randomly generate a pair of point sources Ti and T 2 that 
satisfy a separation condition A > 1/ (2M), which is in 
fact a little smaller than the theoretical constraint, with the 
coefficients of the point sources i.i.d. drawn from the complex 
standard Gaussian distribution. For each Monte Carlo trial, 
we then randomly generate the point spread functions in 
the frequency domain with i.i.d. entries drawn uniformly from 
the complex unit circle, and perform the algorithm by solving 
© using CVX ini- The algorithm is considered successful 
when the Normalized Mean Squared Error (NMSE) satisfies 
LLi pi “ II2 / P*ll2 ^ Pig- E shows the success 
rate of the proposed algorithm over 20 Monte Carlo trials for 
each cell, when M = 8 in (a) and M = 16 in (b), respectively. 
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